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Abstract. We present a theoretical model for the evolution of mixture concentrations in a micro- 
pervaporation device, similar to those recently presented experimentally. The described device makes use 
' of the pervaporation of water through a thin PDMS membrane to build up a solute concentration profile 

inside a long microfluidic channel. We simplify the evolution of this profile in binary mixtures to a one- 

■ dimensional model which comprises two concentration-dependent coefficients. The model then provides 
a link between directly accessible experimental observations, such as the widths of dense phases or their 

, growth velocity, and the underlying chemical potentials and phenomenological coefflcients. It shall thus be 

useful for quantifying the thermodynamic and dynamic properties of dilute and dense binary mixtures. 

PACS. 47.61.-k Micro- and nano- scale flow phenomena - 47.61.Jd Multiphase flows - 64.75.-g Phase 
equilibria - 82.60.Lf Thermodynamics of solutions - 05.70.Ln Nonequilibrium and irreversible thermody- 
' namics 

' 1 Introduction 

] The thermodynamic and kinetic behavior of dense solu- 
' tions is of central interest in a number of industrial and 

■ scientific activities. The so-called "formulation" of multi- 
,: component systems for the production of cosmetics, paints 

and also comestible goods are industrial examples [T] . Sci- 

, entific interests range from mobilities of macromolecules 

. in cells to the statistical physics and material properties 

' of colloidal suspensions [H[3] and other mixtures. In all 

, these applications it is important to know the diffusivity 

' or mobility of solutes, their permeability with respect to a 

I solvent or other solutes and similar non-equilibrium quan- 
tities. Eventually, dense mixtures are complex substances 
which are likely to undergo a number of thermodynamic 
and dynamic transitions of state, which makes their anal- 
ysis a challenge. 

The thermodynamic and dynamic properties of dilute 
and dense mixtures are usually measured by experimen- 
tal techniques such as sedimentation/centrifugation, ul- 
tra-filtration or reverse osmosis [HEIIB] . Depending on the 
properties of the solution under investigation, one or the 
other technique proves to be advantageous. Recently, a 
promising new microfluidic tool for the analysis of aque- 
ous solutions has been introduced, the microevaporator [7l 
[8]. The employed method is applicable to a wide range 
of different solutes, such as electrolytes, surfactants, col- 
loids, or polymers. It makes use of the pervaporation of 
the solvent water through a thin membrane (short vertical 
arrows in Fig.[T^). The solutes cannot pass the membrane 
and build up a concentration profile in an underlying mi- 



crofluidic channel (black dots in Fig.[T^). This proflle of so- 
lute concentration contains information on the transport 
coefficients of the mixture. In experiments, the method 
has been successfully employed for the semi-quantitative 
screening of equilibrium phase diagrams and the determi- 
nation of transport coefficients near equilibrium [5]. As a 
tool for determining phase diagrams, the microevaporator 




Fig. 1. (a) Schematic side view of the micro-evaporator de- 
vice. It comprises a thin PDMS membrane (thickness efa20/im) 
which allows water to pervaporate out of the lower channel 
which is filled by a solution (/i~100/im, L~lcm). (b) More 
details of the flow pattern: The pervaporating solvent induces 
also a concentration profile of trapped solute at the end of 
the channel, here in several dense phases as indicated by the 
different patterns. 
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resembles the functionality of the "phase chip" [IOl[TT|. 
However, the microevaporator allows to go one step fur- 
ther towards the controlled analysis of out-of-equilibrium 
situations. 

The quantitative interpretation of the experimental 
observations in the microevaporator, that is the extraction 
of the transport coefficients, requires a theoretical model. 
The observed quantities are spatio-temporal profiles of so- 
lute densities, the texture of occurring phases, and the 
growth dynamics of dense phases which present a mov- 
ing phase interface. The link between these observables 
on one side and the underlying thermodynamic equation 
of state, the transport coefficients of the mixture, and pos- 
sibly more detailed kinetic parameters on the other side 
can only be achieved by a sufficiently detailed theoretical 
model. At the moment, no theoretical model exists for the 
concentration process caused by pervaporation. We here 
start to fill this gap. 

A model for the phase behavior in microevaporation 
is subject to the following constraints: On one hand, the 
microevaporator device can be applied to very different 
types of solutes, such as colloids, surfactants and elec- 
trolytes. Consequently, a sufficiently general theoretical 
modeling is required which includes phase transitions and 
out-of-equilibrium phenomena namely metastable phases 
and precursor phases. It must further take into account 
the influences of the driving conditions and of the de- 
sign parameters of the channels. On the other hand, the 
model should turn the microevaporator into a quantitative 
tool, making the extraction of thermodynamic and kinetic 
properties of mixtures possible in a systematic and unam- 
biguous way. This aim requires a sufficiently simple model 
which permits to approach several limits analytically, and 
which lends itself to a fast numerical solution of the equa- 
tions. 

The model proposed here can thus be seen as a first 
compromise between the complexity of the subject, includ- 
ing out-of-equilibrium properties, and the requirement to 
allow a fast and direct comparison of its solutions with 
experimental observations. 

In order not to overload the description, we restrict 
it to binary mixtures (water plus a single solute) which 
are described by two spatio-temporal fields, namely the 
concentration of solute 4>{x, t) and a velocity v{x, t) of the 
mixture. The evolution of these two variables is governed 
by equations of the convection-diffusion type, which com- 
prise two transport coefficients, namely D{(j)) for inter-dif- 
fusion of solute and water, and q{4>) for the pervaporation. 
Both coefficients are state-dependent, which is indicated 
by their dependence on the concentration profile 0. 

The structure of the paper is the following: In the rest 
of the current section we give a rough and qualitative de- 
scription of the physics at work, exemplified for a typical 
pattern of dense phases such as the one in Fig. [1)3. This 
specific example will be taken up again in the last section 
of the paper. Section [2] is dedicated to the detailed deriva- 
tion of the model proposed in this paper. It addresses 
in two distinct subsections the modeling of the mixture 
behavior in the microfluidic channel and the pervapora- 



tion through the membrane. All steps and assumptions 
made during the derivation of the final equations are ex- 
plained in the subsections. After the derivation of the gen- 
eral equations we continue to explore several limit cases. 
The section ends with a summary of the model equations 
and the underlying assumptions. In Sec. [3] we focus on 
situations without phase change, where we will study the 
smooth concentration profile of solute in the dilute and 
the dense limits. In the following Sec. [4] we then treat in- 
terfaces between different phases occurring in the mixture. 
The motion of a single moving interface is described first, 
with the aim to extract its dependence on the properties 
of the mixture. Then, we return to the introductory ex- 
ample and analyze the thicknesses of several phase slabs. 
This last analysis is done in a quasi-stationary approxima- 
tion. Both sections [3] and H] employ numerically obtained 
solutions of the model equations and their analytical ap- 
proximations. We finish with a summary and an outlook 
in Sec. [5] 

1.1 The experimental parameters 

The experimental technique and the production of the 
microevaporator device are detailed in Refs. [iMHllQ]. In 
the theoretical treatment in the present paper, we refer 
only to the most essential properties of the device. They 
are displayed in Fig. [TJ We use three geometrical parame- 
ters, the length L of the PDMS membrane, its thickness e 
and the height h of the channel containing the solution. 
The mixture enters the channel from a reservoir with a 
given solute concentration (j)^ ■ This concentration is fixed 
from outside by the experimenter. 

The mixture is driven out of equilibrium by the appli- 
cation of an air flow above the PDMS membrane (thick 
arrows in Fig. [T^) . The air flow guarantees that the perva- 
porated water is constantly removed from the membrane 
and that the permeation continues until the solute "dries" 
out. The out-of-equilibrium situation is thus limited by 
the slow pervaporation. The experimenter's second con- 
trol parameter is the humidity in the air flow. In the case 
of nearly saturated wet air no pervaporation takes place, 
because the chemical potentials of water on both sides 
of the PDMS membrane will become equal. In general, 
the two external "knobs," namely the reservoir concentra- 
tion (pL and the air humidity need not be constant but 
can be changed according to a temporal protocol. 

Due to the slow pervaporation process, a typical ex- 
periment takes hours to days. According to Ref. [7], the 
intrinsic time scales of the device are in the range of sec- 
onds (for the evaporation and the mechanical adjustment 
of the membrane) up to fifteen minutes (for the diffusion) . 
We may therefore hope that the mixture is not too "far 
from equilibrium," which is a necessary requirement for 
the following description. 

During the evaporation, several quantities can be mea- 
sured, such as the total amount of mixture which passes 
through the device during the experiment. It gives some 
information on the amount of solute in the channel and 
on the rate of water pervaporation. Within the channel. 
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more information on the mixture can be gathered: If there 
are several phases in the mixture and if they can be dif- 
ferentiated e.g. by their texture, then the position of the 
interfaces can be tracked as a function of time. This has 
been done for the case of the surfactant AOT (docusate 
sodium salt) in water in Ref. [8] and for a salt solution in 
Ref. [7]. More information is obtained if the full density 
profile can be measured as a function of time 7j- These 
types of observables are in the focus of the present theo- 
retical description. 

1.2 An intuitive example 

The central interest of the present theoretical description 
is the precise functional forms of two coefficients, one of 
which is the inter-diffusion coefficient £'(0) of solute and 
solvent, the other one q{(j)) quantifies the pervaporation of 
solvent through the PDMS membrane. Both coefficients 
do depend on the local concentration (j>{x, t) of solute. All 
thermodynamic and kinetic properties of the mixture are 
implicitly contained in these two coefficients. In particular, 
they can be expressed in terms of the chemical potentials 
and of the kinetic properties of the mixture, both of which 
depend on 0. 

Suppose the situation in Fig. [TJ), which has been ex- 
plored experimentally using the surfactant AOT (docusate 
sodium salt) in water |8j: The channel is filled by four dif- 
ferent phases of the mixture, which are ordered by their 
solute concentration. Three of them are dense, indicated 
by different patterns in Fig. [T|d, the fourth phase is dilute. 
The pervaporation of water, quantified by introducing the 
coefficient q{4>) gives rise to a flow from the reservoir (long 
curved arrows) which compensates the water loss at the 
PDMS membrane. On its way through the phases the wa- 
ter drags the solute with it. It thereby induces the concen- 
tration of solute, and, at higher and higher concentration, 
the appearance of the three dense phases which would be 
absent without flow. 

The final concentration profile is thus a dynamic equi- 
librium between convective and diffusive processes: Within 
each phase slab the concentration gradient gives rise to 
diffusion which in turn tries to shrink the phases. While 
the diffusion emerges together with the concentration gra- 
dients, the convective transport, or "drag" is reduced by 
the solute becoming dense. The solvent has to permeate 
through the dense phases of solute before it reaches the 
membrane and pervaporates. The denser the solute, the 
less solvent may pass. Furthermore, the more the water 
flow tries to compress, or to "squeeze" the dense phases, 
the harder it becomes to squeeze them even further. One 
or the other of the mentioned mechanisms may remind the 
reader of terms such as "osmotic compressibility," "diffu- 
sion," "mobility," "permeability," depending on the liking 
and the experience of the reader. In fact, all these effects 
are linked to one another, and we will provide the connec- 
tion between some of these terms in Sec. 12.31 We thereby 
also show the equivalence of these terms, in so far as they 
can expressed one by the other. In our model we abandon 
the use of terms such as mobility, osmotic compressibility 



and permeation in favor of the diffusion coefficient D((j)) 
and the coefficient of pervaporation q{(j)), which are ex- 
pressed in terms of the chemical potential of water and of 
the phenomenological coefficient for inter-diffusion. 

An important observation in the situation of Fig. [T]d 
is the following: As more and more solute is transported 
from the reservoir, all three dense phases grow. It has 
been observed experimentally that the two interior "sand- 
wiched" phases apparently grow at constant thickness. 
They loose as much solute as they gain, while the last, 
the most dense phase continues to grow. Keeping in mind 
the dynamic equilibrium between convective and diffusive 
processes within each phase we therefore find in the thick- 
nesses of these middle phases a fingerprint of the diffusion 
coefficient and of the pervaporation coefficient. 

One may now ask whether a large diffusion coefficient 
will finally lead to a thinner or to a thicker slab, as com- 
pared to one with a smaller diffusion coefficient? On one 
hand, the diffusion coefficient scales the gradients of the 
concentration profile which try to dissolve the phases. It 
should thus lead to small thicknesses of the quenched 
dense phases. On the other hand, for fixed values of the 
concentration at the phase interfaces, a large diffusion co- 
efficient will render the profiles fiat and therefore lead to 
thick phases. This may appear paradoxical on a first sight, 
but it only reflects the same feedback mechanism as be- 
tween flow of solvent and permeability. The resolution is 
that we really need the full concentration proflle and the 
velocity fleld in order to understand the phase thicknesses. 

2 Detailed derivation of the model 

In this section we derive a set of equations which model 
the evolution of a binary system, reduced to one spatial 
dimension x along the evaporation channel. As variables of 
state we use a concentration (/)(a:, t) of solute and a mixture 
velocity v'^{x, t). Special focus will be put on their precise 
meaning in the evolution equation that will be derived in 
the following sections, 

= -i'»g(0), (1) 

dtc^ = -d^[<l>v^ -D{c^)d^c^]. (2) 

Here, denotes the specific mass of water; dt and dx are 
partial derivatives with respect to time and space. The 
physical content of the equations is the independent bal- 
ance of solvent and solute mass. The variables are chosen 
such that q{(j)) reflects the loss of solvent (water) via per- 
vaporation, and D[(j)) is a state-dependent inter-diffusion 
coefficient, which contains information on the thermody- 
namic and dynamic properties of the mixture. The expres- 
sions for D(0) and g(0) will be provided by equations ([3T|) 
and ((46|) below. 

2.1 Diffusion model for the binary solution 

The aim here is to state clearly all underlying assumptions 
which lead to the model equations ^ and ([2]), starting 
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from mass conservation. This procedure will allow to point 
out the restrictions of the model. We employ the theory of 
linear out- of- equilibrium thermodynamics as summarized 
in the book by de Groot and Mazur [12] . We try to adopt 
the notation used there, also for an easier comparison with 
the results by Peppin etal. [4] for the link to sedimenta- 
tion, ultra-filtration, and reverse osmosis. 



2.1.1 Balance equations for binary mixtures 

The balance equations for the mass densities ps(r, i) and 
pu,(r,t) of solute (index s) and solvent (index w) read 



dtPk = - div(jfc), with k e {s, w}. 



(3) 



The mass current densities jfc are split into an advected 
and a diffusive part. 



(4) 
(5) 



The Ji denote the "diffusive currents" , simply defined by 
the difference pi(vi — v). A common choice for the mixture 
velocity is the mass-average of the individual constituent 
velocities v^, weighted by the mass fractions c^, 



with 



— and 



k 



(6) 
(7) 



This choice makes the two diffusive currents opposite to 
each other, 

Jw = Js- (8) 

Formally, the mass densities are taken to be functions 
of the local thermodynamic variables of state as Pi{T, P,Cs). 
The use of local temperature and pressure fields shows 
that here and in the following we will assume local ther- 
modynamic equilibrium. The evolution equations ([3]) are 
accompanied by the following constraint, which is an im- 
mediate consequence of the chemical potentials being in- 
tensive quantities (or of the Euler relation). 



1 = t^sPs 



VwPu 



(9) 



The specific volumes per mass, Vi, are defined as the deriva- 
tives of the chemical potentials per mass, 7i(r, P, Cg), with 
respect to the pressure. 



u,{T,P,Cs) :-- 



l^{T,P,Cs) 



Pi{T,P,Cs) 



(10) 



(11) 



Here, pi is the chemical potential per molecule, and 
the mass of such a molecule. The derivative in Eq. (|10p 
equals the change of volume induced by adding an in- 
finitesimal amount of mass of type i upon keeping temper- 
ature, pressure and all remaining masses constant. Equa- 
tion ^ can be used to introduce the volume fractions 



(j)i{T, P, Cs) := ViPi- Below, equation ^ will also be the 
point of entrance for the notion of simple mixtures and of 
incompressibility of the constituents. Apart from Eq. ([9]), 
another consequence of the intensivity of chemical poten- 
tials is the following relation between the derivatives of the 
two chemical potentials with respect to the mass fraction 
of solvent, c,, 



T,P 



dc. 



(12) 



T,P 



This equation establishes the standard logarithmic behav- 
ior of the chemical potentials in the limits Cg ^ and 

Together, equations ^ and (fT2|) reduce the number of 
relevant chemical potentials to one, as one can always be 
expressed by the other. This fact is important for the num- 
ber of kinetic coefficients within the theory of linear non- 
equilibrium thermodynamics. Applying this framework in 
the present description, the diffusive current is taken 
to be proportional to the gradients of chemical potentials. 
The proportionality factor is called the phenomenological 
coefficient C. Using the relation of Gibbs-Duhem, only the 
difference of chemical potentials can occur, and we receive 
for the diffusive current 



-C V(7s - 7t, 



(13) 



We do not present the details of the derivation of this 
equation here, as it implies a rather tedious consideration 
of all possible symmetries in the entropy production. In- 
stead, we refer to Eq. (IV. 15) in the book by de Groot 
and Mazur |12j or to the summary in the paper by Pep- 
pin et al. j4] . Note that external forces such as gravity can- 
cel out and do not appear in equation ([T3|). This fact is 
due to the chemical potentials being defined as energy per 
mass rather than per number of particles. Furthermore, we 
have omitted a term proportional to the temperature gra- 
dient, as from now on we take the temperature as uniform 
and constant. We will omit it also from the argument lists 
of thermodynamic functions. The phenomenological coef- 
ficient C may in principle depend on all thermodynamic 
variables of state as well as on their derivatives, and even 
on the mixture velocity v and its derivatives. However, 
if the assumption of local thermodynamic equilibrium is 
well satisfied, the dependence on the mixture velocity field 
should be negligible, such that we use a function C{P, Cg) 
for the moment. The chemical potentials are also func- 
tions of the two thermodynamic parameters: 7s (P, c^) and 
7u,(P, Cs). This functional dependence allows to separate 
the influence of pressure gradients from concentration gra- 
dients in Eq. P^ . 



Js = -C 



{vs-Vw)VP- 



d{ls - Iw) 



dcs 



T,P 



(14) 



With the aid of relation the diffusive second term can 
be rephrased using either of the chemical potentials. 



dcs 



1 d-is 



T,P 



1 - Cc dCa 



1 d^u 



T.P 



dcs 



T.P 



(15) 
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Equation (fT4| has been used by Peppin etal. 4, Eq. (28)] 
to demonstrate that Fick's and Darcy's laws are limiting 
cases of the same equation. Below, we will use it to de- 
fine the diffusion coefficient and to make the connection 
between the different methods micro-pervaporation, sedi- 
mentation, and ultra-filtration. 



2.1.2 Reduction to ID 

We shall now disregard the details of the flow pattern in 
the channel and reduce the three-dimensional evolution 
equations ([3]) to one dimension. We do this by averaging 
over the two spatial dimensions perpendicular to the chan- 
nel axis, as depicted in Fig.[2l For a given uniform height h 
and width b of the channel we obtain from integration over 
a slice SI with area |J7| = bh the one-dimensional counter- 
parts of ([3]), 



dtPs 
dtPw 



-dxjs 
-dxjv 



The new densities and currents turn out to be 



Pk{x,t) := 
jk{x,t) := 



1 

W\ 



dy dz pk{x,y,z,t), 
dydz • jk{x,y, z,t). 



(16) 
(17) 



(18) 
(19) 



The water loss due to pervaporation through the PDMS 
membrane is all included in q{x,t). It represents the inte- 
gral of water current in normal direction over the bound- 
ary of the channel. 



q{x,t) :-- 



1 

W\ 



N • j^{x,y,z,t) 



(20) 



b 



i-^ N-jUx,y.Kt)dy. (21) 



The pervaporation has a dominant contribution from the 
upper boundary where the thin layer of PDMS allows a 
steady water pervaporation. In addition to equations (jl6p 
and pT|) we expect an averaged version of constraint ([9]) 
to hold. In the same manner as above in Eqs. ^ and ^ 
we may subdivide the total currents into advectivc and 
diffusive parts. 



2.1.3 Incompressibility and simplicity of mixtures 

We now assume both constituents of the mixture to be 
individually incompressible, having constant specific vol- 
umes Vi. This assumption will allow to reduce the set of 
three equations p6)) . p?)) . and ([9]) to only two equations. 
The assumption, however, is quite strong as it implies the 
mixture to be both incompressible and simple. 



dl/p 



dP 



dl/p 



dcs 



= const. (22) 




t t t t<li^,t) 



^= ja{x,t),jn,{x,t) 

Fig. 2. Visualization of the averaging used for reducing the 
balance equations to one dimension, (a) The area of integra- 
tion Q in three dimensions, (b) The one-dimensional represen- 
tation. 



The term simple mixture here means that the total mass 
density is a linear function of the mass fraction. It implies 
linear relations between mass density, mass fraction and 
volume fraction. We will therefore use the terms "density" , 
"fraction" and "concentration" interchangeably. Simple 
mixtures are for example glycerol in water }13j or rigid 
spheres in an incompressible solvent of much smaller par- 
ticles. The constraint ([9]) allows to write one density as a 
function of the other, — p^ips)- The evolution equa- 
tion P?|) for Pw{ps) then becomes 

{pw{Ps) - ^;^P^dxV + dxJw - ^p^dxJs = -q- (23) 
V dps / ops 

This equation can be understood as a differential equa- 
tion for the velocity field, similar to the one we are look- 
ing for llj. However, the occurrence of the diffusive cur- 
rents Js and makes it unfavorably dependent on the 
concentration variable, see equations (fT4|) and ([8]) . A more 
convenient equation is obtained by using from the begin- 
ning the volume-averaged mixture velocity 



(24) 



instead of the mass-averaged one. The total currents are 
then written with modified diffusive currents, 



Jk 



^Pkv° + Jk, ke{s,w}. 



(25) 



With the volume-averaged mixture velocity, and by set- 
ting the specific volumes constant, one obtains as relation 
between the new diffusive currents. 



J?., - J° (26) 



instead of equation ([5]). This property makes the diffusive 
currents vanish from the differential equation for u°(x,t), 
which now comes close to the form we seek ([IJ , 



dxv° = -Vwq{P,(t>)- 



(27) 



In terms of the volume-averaged mixture velocity, and 
with constant v^, the volume fraction of solute evolves as 



T,P 



dt(t>s = -dx <PsV° + {vw<Ps + Vs{l~(l)s)) Js 



(28) 
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Or, omitting the index s and collecting also equations 
and ()15p . we finally arrive at 



bv^ ~K{P,4,)d^P-D{P,<f>)d^ 



(29) 



with the shortcuts 



D{P,(I} 



[i^s + 4>{v^n - ys)\ C{P, 0) {us ~ yw), (30) 



(31) 



Equations and with the coefficients q and D 
from (|^D1) and (|3ip come close to the model we announced 
in the introduction. They nicely separate the effect of per- 
vaporation, which is all contained in the coefficient q, from 
the bulk thermodynamic and kinetic properties of the mix- 
ture. These properties are implicitly contained in the func- 
tional dependence of C and 7^ on P and on (f>. 



2.1.4 The role of the pressure gradient 

With equations (P7|) and (^5)1 we have found two equa- 
tions for yet three variables v'^, 4>, and P. We thus have 
to seek for a closure of the above set of equations. The 
contribution by the pressure requires some further atten- 
tion. The pressure occurs as an argument in several func- 
tions, but also its gradient adds a contribution to the dif- 
fusive current in Eq. This gradient appears also in 
the (compressible) Stokes equation. In principle, it would 
have been possible to assume local mechanical equilibrium 
already in the 3D theory, such as the a:-component of the 
compressible Stokes equation. 



= -d^P + T]Av^ + (A + T]/3)d.j, divv 



(32) 



and then reduce it to one dimension together with fit- 
ting boundary conditions. The shear and bulk viscosities 
are denoted by 77 and A, respectively. However, this proce- 
dure raises fundamental difficulties: We showed above that 
the natural description of pervaporation of incompressible 
mixtures involves the volume-averaged velocity v^. The 
natural formulation of the Stokes equation uses the mass- 
averaged velocity. Only for the special case Vs = i^w both 
velocities do coincide — but in this very case the term with 
dxP vanishes due to X = 0. When expressing v by v°, the 
term with dxP in Eq. ([29| will contribute to both other 
terms in the expression for the diffusive current in (j29p . 
We do even expect a coupling of velocity and density vari- 
ables, as well as a second derivative (9^0 due to the second 
derivatives in the Stokes equation ([35)1 . 

The precise analysis of this splitting is beyond the 
scope of the present paper. Here, we have to restrict the 
validity of the model to cases where the influence of the 
pressure gradient is negligible. An estimation |14j indicates 
that this restriction is indeed weak: In the dilute limit, 
using values from the experiment |7j, the term with the 
pressure gradient in Eq. ([29]) is several magnitudes smaller 
than the diffusion term with the concentration gradient. 
In what follows we therefore use the simplistic model with 



a homogeneous pressure variable (which we omit in the 
notation from now on): 



dxV° = -iywq{<f>), 



(33) 
(34) 



with the evaporation coefficient q{(j)) and the diffusion co- 
efficient D((j)), both depending only on the volume frac- 
tion (j). 

As we take the microevaporation channels to be hor- 
izontal, our assumption of setting the pressure gradient 
to zero coincides with the assumptions employed by Pep- 
pin et al. 4 and also in the context of sedimentation [5] 
where the pressure gradient balances gravitational forces. 
It therefore opens the possibility to consistently compare 
with their model equations, see Sec. 



2.1.5 Boundary conditions 



Equations f33|) and (|34p are to be accompanied by an 
appropriate number of boundary conditions. For the hy- 
perbolic differential equation p3p we require one, for the 
parabolic equation (j34p we need two boundary conditions. 
At the end of the channel, both the mixture velocity and 
the total solute current vanish. 



^^°(0,t) = 0, 
9,(/.(0,i) = 0. 



(35) 
(36) 



As a third boundary condition we use a fixed concentra- 
tion at the inlet of the channel, reflecting the concentra- 
tion in the reservoir. 



<t>{L,t)^(l>i 



(37) 



The model equations could equally handle a time-depen- 
dent reservoir concentration 4>L{t)- In a setup with a fixed 
total amount of solute in the channel, it is mathematically 
also possible to replace the latter boundary condition by 

an integral condition setting 4>{x, t) dx to a given value. 



2.1.6 Moving interfaces between phases 

In the introductory example we mentioned that the po- 
sitions of phase interfaces and their velocities are conve- 
nient observables in the pervaporation experiments. Such 
experimental observations have been reported in Ref. [7] 
and in more detail in Ref. [8] for a solution of the surfac- 
tant AOT (docusate sodium salt) in water. The essential 
observation in the latter reference is the decrease of the 
interface velocity with time, which will be investigated in 
more detail below in section W7\\ A quantitative analysis 
of this decrease can give insight into the mechanisms at 
work, especially on the relative importance of diffusion 
and advection caused by pervaporation. 

Throughout this paper, we model the interfaces within 
the framework of local equilibrium thermodynamics, which 
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X{t) X 

Fig. 3. Sketch of the modeling of sharp interfaces between 
different phases. The concentration values left and right of the 
interface are fixed by a Maxwell construction on the chemical 
potential (upper panel), the two regions shaded in gray have 
equal area. 

implies that the interfaces are infinitely sharp and that for 
each such interface there are two unique concentration val- 
ues at its left-hand and its right-hand sides. These values 
are determined by a Maxwell construction on the local 
chemical potential function. Figure [3] sketches this con- 
struction and the resulting concentration profile around 
the interface. The fixed concentration variables serve as 
two Dirichlet boundary conditions for the evolution equa- 
tion 

0(X(^),t) 0(xW,t) -(/.W. (38) 

The superscripts [t) and (r) indicate the left and right 
limits towards the interface position X{t). 

The conservation of solute mass, which is valid also 
across the interface, provides us with a continuity condi- 
tion at the interface. If the boundary does not move, the 
solute currents at the left-hand and at the right-hand sides 
of the interface must coincide. Any possible difference be- 
tween these incoming and outgoing currents must lead to 
a growth of one phase and thus to the movement of the in- 
terface position X{t). Its velocity times the concentration 
jump at the surface must equal the difference between the 
currents, 

■ ^^^^ 

We further require continuity conditions for the veloc- 
ity field. As we did in the derivation of the bulk equations, 
we also here want the mixture velocity to reflect only the 
influence of the pervaporation. As we do not expect the 



pervaporation to diverge near phase interfaces, the veloc- 
ity field should be continuous across such an interface, 

v°{X^^\t)^v''{X'^'-\t). (40) 

This continuity equation is consistent with the conserva- 
tion of solvent mass across the interface. 



2.2 Pervaporation through the PDMS membrane 

The pervaporation of water through the PDMS layer can 
be described on the same level as the advective and diffu- 
sive currents in the evaporation channel. We now focus on 
the binary mixture of water and PDMS, as solutes cannot 
enter the PDMS domain. The notation is such that the 
properties inside the PDMS membrane will be denoted by 
an overbar. As we did above, we also here split the total 
current density ]w of water into convective and diffusive 
terms. This time it is more convenient to take the velocity 
of PDMS, which is zero, as the reference velocity of the 
mixture. The mass current density of water in the layer 
then becomes ^ 

]w = . _ - Jtt,, (41) 

where c^, is the mass fraction of water inside the layer, 
and Jk, is the corresponding diffusive current. In com- 
plete analogy to equation (|13p . the diffusive current can 
be written in terms of a phenomenological coefficient C, 

I. = -^^^^V(7.-7pd). (42) 

with the chemical potential of PDMS denoted by 7pd . We 
understand the pervaporation as a diffusive process and 
again omit the pressure-dependence of the chemical po- 
tentials. Then, the water current is proportional to the 
gradient of water concentration. The diffusion process is 
further assumed to be sufficiently established that the wa- 
ter current does not vary within the membrane. In the 
dilute limit, where the chemical potential exhibits a loga- 
rithmic behavior, 

kT 

7» (T, P, ) = 7: (r, F) + In , (43) 

the uniform water current then corresponds to a constant 
gradient of water concentration inside the membrane. The 
water current through the membrane is then proportional 
to the concentration difference across the membrane of 
thickness e, 

_ T _ c^,{z=h+e) - Cw{z=h) C{P,Cw) dj^i 

e E„^o (1 - Cu,)2 dcw 

(44) 

To make the connection with the conditions in the ad- 
jacent channels of the mixture and of the air ffow, we 
assume that the chemical potentials at both sides of both 
boundaries of the PDMS membrane have matching values, 

7t„(z=/i) = 7t„ and j^{z=h+e) ^ . (45) 
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The evaporation coefficient q can now be expressed as a 
function of the chemical potential jwi^) in the mixture 
channel. From equations (|44)) . ([43]) . and its definition (|2ip 
we finally obtain the following dependence of the pervap- 
oration coefficient on the volume fraction, 



with the diffusion coefficient 



1 D 

v,,, eh 



(m^ airl 

exp|^7.W|-exp|— 7. I 



(46) 



We have collected several of the above constants in the 
shortcut Z), 



n ( rn^-^ \ y ^ 97„ 



(47) 



Expression (j46p contains the further assumption that the 
chemical potential directly at the interface between mix- 
ture channel and PDMS layer is the same as the average 
chemical potential in the mixture which we obtained from 
the reduction to one spatial dimension. In the present 
framework of a one-dimensional description, there is no 
alternative to this approximation. 

At the present point we have completed the deriva- 
tion of the model, consisting of the two differential equa- 
tions ([55]) and for the volume fraction (p and the vol- 
ume-averaged velocity field together with the concen- 
tration-dependent diffusion coefficient (|3ip and the evap- 
oration coefRcient pS)) . 

2.3 Connection with sedimentation and ultra-filtration 

We now make the connection between the two coefficients 
D{(f)) and q((/)) and other pairs of coefficients. This will al- 
low to compare different techniques of measuring chemical 
potentials and phenomenological coefficients and to verify 
the consistency of the outcomes. We will not make use of 
these connections in the present paper. The equations in 
this section can rather be seen as a convenient reference 
for the reader working on different techniques. However, 
the comparisons give some more intuition on the variables 
and coefficients which are used. Above, in equations (|3ip 
and (|46|) we have already found the link to the pair 7(0) 
and £((/)). We continue by expressing also the coefficients 
of sedimentation and of permeability in terms of 7('/>) and 
C{(j)). Each of these two coefficients forms together with 
the diffusion coefficient a pair which allows to determine 
7(0) and C{(j)) and thus to make the fink to our original 
coefficients D{(f)) and q{4)). 



2.3.1 Sedimentation factor and osmotic compressibility 
factor 

In the context of sedimentation of colloids, Russel, Sav- 
ille, and Schowalter express the thermodynamic and 
dynamic properties in terms of the sedimentation factor 
K^^^{(j)) and the osmotic compressibility factor Z^^^{(f>). 
Both can be introduced by rewriting Eq. as follows [31 
compare Eq. 12.5.6], 

dtct) ^-div(^<t)v" +Vo4>K^^^{4>) ~ D{(t>)V^^ (48) 



D{<j>)^D,K^^^{^)^[cpZ^^H4)] (49) 



and with the shortcuts 
kT 



Un 



2a 



1 1 

977 Vj/„ r 



(50) 



with r] denoting the viscosity of pure solvent, a the radius 
of a colloidal sphere, and g the acceleration by gravity. The 
volume-averaged velocity in Eq. (|48p vanishes by construc- 
tion of the sedimentation experiment in a closed container. 
In the comparison with Eq. (|29p the remaining two terms 
thus lead us to the sedimentation coefficient i^^ss gjygjj^ 



pg 2a^ V 



-I- 



(51) 



The osmotic compressibility factor, originally defined by 
the osmotic pressure of solute can equivalently be 

expressed in terms of the chemical potential of water, 

Z^^^(6) = ^ "^^^ = 3 7^-(0) -7^-(0) /.2) 
47ra3 #T An u^(j>kT ' ^ ' 

Peppin, Elliott, and Worster [4] define the sedimenta- 
tion coefficient S^^^ in a marginally different way. Their 
definition leads to the following variant of the transport 
equation ^ (see Eq. (32) of Ref. gj), 

dt(t) = - div((/)v" + g4>K^^'^{(P) - D{<t))V(p^ (53) 

In their treatment, the pressure gradient is assumed to 
balance the force density by gravitation, VP — pg. In our 
notation of Eq. ([50)) , their sedimentation coefficient reads 



S 



PEW 



1 



(54) 



2.3.2 Permeability 

For the coefficient of permeability k^^^ Peppin et al. give 
the resuh (Eq. (37) of Ref. g]) 



^PEW 



(55) 



Except for the dilute limit, the phenomenological coef- 
ficient £((/>) can be understood as a direct measure of 
the permeability of the solute with respect to the solvent. 
This will be useful below when choosing a model function 
for £(0). In particular, Eq. ([55]) helps choosing the value 
of £(1), which has influence on the qualitative behavior 
of the solutions, see the discussion of Fig. [51 In the oppo- 
site limit, Eq. ([5T|) indicates that £(</>) oc for small 
This last requirement leads to a finite non-zero diffusion 
coefficient in the dilute limit, see below in Fig. [T3l 
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2.4 Summary of the model equations, and their 
preparation for a numerical implementation 

We here summarize the model equations ([1]), ([5]), ([51]) . 
and (|46)) for their further use in the following sections. We 
repeat their boundary conditions and the condition for 
the interface movement. We also reformulate the model 
equations in such a way that they may be discretized and 
implemented more easily. The dimensionless equations of 
evolution ([T]) and ^ read 



'dx'j' = -dx' 



He 



(56) 
(57) 



They describe the evolution of an incompressible binary 
mixture in one effective spatial dimension under the fur- 
ther assumption of homogeneous temperature and pres- 
sure. The variables carrying a prime are made dimension- 
less by the following scales for length (x), time (t) and 
the mixture velocity (w°). The volume fraction (p needs no 
further rescaling. The scales have been chosen such that 
all dimensionless variables are of the order unity as far as 
possible, 

x = Lx', t = t' /{v^q{Q)), v° ^ Lu^q{0) v°' , (58) 
g(0) = g(O)g'(0), Di<P) = DiO)D'{^), (59) 



Pe 



(60) 



Here, the Peclet number Pe is a global Peclet number of 
the system. The scales for the diffusion coefficient and the 
coefficient of evaporation, which are simply the coefficients 
evaluated at zero concentration density (£'(0) and q{Q)) 
are chosen such that Pe presents an upper bound for the 
true (local) Peclet number Peioc{x,t) which can only be 
determined a posteriori from the solution. In numerical 
solutions we have found the local Peclet number to be 
always maximal at the entrance of the channel and unity 
at the end. However, it does not need to be monotonous. 

The dimensionless version q'{4>) of the pervaporation 
coefficient is, according to the expression we found in ()46p . 



exp 



/ w 



exp 



'{^(7^0) 



It depends on the volume fraction only via the chemical 
potential of water 7^,(0), all other parameters are con- 
stants or external driving parameters. The dimensionless 
diffusion coefficient reads, according to equation pTjl . 



D'{4>) - 



907m, (0) lim C{ip)/ip' 



(62) 



It obtains it dependence on </> from ^^{4') ^^J^d from the 
phenomenological coefficient £(0) for inter-diffusion. 

The boundary conditions ([55)1 . ([5^ . ([57)1 read in their 
rescaled version, 



,,0/ 



(x',t') = ata;' = 0. 



(63) 



dx,(l){x' ,t') atx'^O, 
4>{x ,t') — (f>L at x' = 1, 



(64) 
(65) 



as well as the boundary/continuity conditions at moving 
phase interfaces (pS)) and (|40|) . The velocity ((39|) of an 
interface becomes 



Pe [cj)W - (l){r)] 



(66) 



The following sections will require solutions for the 
model equations, partly for arbitrary functions q{(j)) and 
D((f)). Such solutions are available only by numerical meth- 
ods. For the numerical discretization we transform equa- 
tions ((56|) and (|57p into time-dependent systems of co- 
ordinates (,{x',t') which range from to 1 in each phase. 
x'(^, t') interpolates linearly between the (moving) bound- 
aries, 

x\^,t') = X"^'\t') + ^[X'^^\t') - X"^''\t')] (67) 

Here, X'^^^ and X''^^^ are the positions left and right of 
the phase in question. We further use the shortcut AX' :— 
X'('^'> —X''-^\ The coordinate transformation adds a second 
advection-like term with velocity 9fx'(^,t') to the evolu- 
tion equation, such that we obtain for the model equations 
in each phase. 



5t'0(e,t') 



dtx'iW) 



AX' 



8(6 - di 



Z\XV(0), 

D'{(f)) 



(68) 



AX' 



Pe{AX'Y^^^ 



(69) 



Equations (pB]). and are then solved using the 

method of lines, with a finite-volume discretization and 
linear interpolation of (f> in space. The time-stepping is 
done by a standard explicit fourth-oder Runge-Kutta al- 
gorithm [TS]. 



(^52) 3 Prototype solutions: a single phase 



In this and in the following sections we discuss several 
different solutions of the above model equations, which 
are ([1]) and ([2]), completed by the coefficient functions g(0) 




Fig. 4. Sketch of the typical single-phase situation described 
in Sec. [3] 
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and D{(t)) from Eqs. ([3T|l and (|46|) — Or, in their dimen- 
sionless versions, Eqs. dM]), dM]), and We here 

choose two functions q{(j)) and D{(j)) and solve the equa- 
tions numerically, accompanied by analytical approxima- 
tions. The aim is to display typical solutions and to ex- 
plore the variety of solutions by using different functions 
q{(j)) and D{(f). These solutions are meant to guide fu- 
ture experiments and their interpretation, where the final 
target is to employ the model for extracting the functions 
q{4>) and -D(0) from experimental data and to invert them 
to obtain the chemical potentials and the phenomenologi- 
cal coefficients. This will then require to solve the inverse 
problem of what we do in the present paper. We start 
with the description of three typical situations with a sin- 
gle phase and smooth concentration profiles, see also Fig.|4] 
for a sketch. 



3.1 Solvent only 



0.1 r r 




Fig. 5. Growth of a dilute solute density according to the two 
regimes "Gaussian" and "hyperbolic ramp." The plot shows 
the numerical solution with parameters (g' = 1, D' = 1, Pe = 
100, solid curves) together with the analytical approximation 
from Eqs. and (I72|l (dashed curves), see also Ref. 0. 



For setups containing only water and no solutes, the water 
loss is the constant (/(O), which has been measured experi- 
mentally by the total consumption of water per time . A 
situation with water only therefore permits to identify the 
constants in expression (|46p for the pervaporation coeffi- 
cient. A corresponding series of experiments with varying 
parameters can be used to verify that the assumptions we 
made for are indeed fulfilled. Possible variations can 
either be geometrical such as changing the length L or 
the height h of the channel, or the membrane thickness e. 
Other variations concern the driving parameters, such as 
the solute concentration 0^ in the reservoir, or the princi- 
pal control parameter during the experiment, which is the 
external chemical potential 7^"^ of water in the air flow. 
This chemical potential is controlled by the humidity of 
the air circulating through the top channel. Equation (|46p 
predicts a linear dependence of the pervaporation on the 
specific humidity c^"', with an offset. 



at Peclet number equal one. Both regimes can be identi- 
fied in Fig. [5] where a numerical solution of equations ^ 
and JH) is depicted together with the analytical approxi- 
mation from Ref. These approximations are given by 
the solution of 



j'(l) = j'{x') « V 



^q'x'q 



/(I) 
q'x' 



(71) 



for the hyperbolic ramp, and by the solution of 

D' 



= fix') = 



-q'x'c, 



Pe 



(72) 



for the Gaussian profile. The width of the Gaussian hump 
thus scales as L/\/'Pe. The values given in Ref. 7 indicate 
a Peclet number of 100. 



ID 

m = — r 

V,,, eti 



where the slope of the line and the offset may in prin- 
ciple depend on the pressure. Quantitative information 
on this dependence would thus offer a validation of the 
current model and yield insight into the thermodynamic 
mechanisms inside the PDMS layer. Previous studies [T^ 
assumed that the pervaporation is driven by a gradient of 
concentration and not of pressure, as we did it in Sec. 12.21 



3.2 Dilute solutions 



3.3 From dilute to dense single-phase solutions 

The mere model equations ([1]) and ([2]) include no mech- 
anism to prevent unphysical values <f>{x,t) > 1. We thus 
have to ask what mechanisms the model provides to en- 
sure (j) < I'm the dense limit when (p approaches unity. For 
dense solutions, both coefficients q{4)) and -0(0) exhibit 
non-trivial dependences on (j) which provide two indepen- 
dent physical mechanisms to keep (j) below unity. These 
mechanisms are described in the following paragraphs. 
Numerical examples for both are depicted in Fig. [51 



The situation with dilute solutions, in which both the 
evaporation coefficient q and the diffusion coefficient D are 
essentially constant, has already been explored in Ref. [7]. 
Two different regimes were identified in the solution, first a 
hyperbolic ramp corresponding to a steady uniform current 
at high Peclet numbers, and second a Gaussian growing 
linearly in time, corresponding to the accumulation zone 



3.3.1 Vanishing pervaporation 

One mechanism for keeping the volume fraction bounded 
is governed by the pervaporation coefficient q{4i). Equa- 
tion (I46p for this coefficient implies an equilibrium value (jf^ 
at which the pervaporation vanishes. If the concentration 
exceeds this value, the pervaporation direction is inverted. 
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which stabiHzes 4> around the equihbrium value, see Fig. [71 
In combination with the boundary condition ()35p for the 
velocity field we find a flattening of the concentration pro- 
file at (j) = (f)'^'^, at least in an interval starting with x = 0. 
In this interval, the velocity remains approximately zero. 
A numerical example of such a profile is given in Fig. ^p. 

3.3.2 Diverging diffusion coefficient 

The other mechanism leading to volume fractions t) 
which strictly cannot exceed unity is that the chemical po- 
tential of water diverges at ^ 1. In this case, the solvent 
becomes dilute in the solute, leading to the standard log- 
arithmic form of the chemical potential of water in the so- 
lute concentration (cf. Eq. (fT2|) or Eq. ([43|l ). Together with 
the chemical potential also its derivative diverges in the 
limit — > 1. In other words, the osmotic compressibility 
of the solute phase diverges as well: The more the solute is 
concentrated, the more difficult it becomes to concentrate 
it even further — and it cannot be concentrated beyond a 
value t/)**^ which is given by thermodynamics as the point 
where the chemical potential diverges. 0**^ will be unity in 
most cases. 

Mathematically speaking, the diffusive currents tend 
to flatten the spatial concentration profile. In the limit 
of an infinite diffusion coefficient — which is proportional 
to the derivative of the chemical potentials — the volume 
fraction (f>(x,t) become arbitrarily flat. It will then not 
grow beyond the value 0''* at which the diffusion coeffi- 
cient diverges. A numerical example for such a situation 
is depicted in Fig. [5^. This flattening poses a bounding 
mechanism which is absent if the diffusion coefficient stays 
finite or even vanishes. 

It depends on the values of (f)'^'^ and 0*'^ which of the 
two mechanisms takes place first in a given mixture. In 
the present description, however, which is based on ther- 
modynamics, the logarithmic divergence must be located 
at f/)'*^ = 1 which is always larger than the equilibrium 
value In a mixture of small solutes, which are of 
molecular size, the value (j)^'^ = 1 should indeed be ap- 
proachable. It is not approachable in a system of rigid 
spheres in a much smaller solvent, where there is a close- 
packing limit 0*'* < 1. In such systems with large solute 
particles, the solvent can still pass through the holes be- 
tween the spheres, and the pervaporation never ceases [T7] . 
We acknowledge that such a system cannot be described 
by our model — a deficiency which stems from the treat- 
ment of the pressure, and which is therefore shared by 
any other thermodynamic model that assumes mechani- 
cal equilibrium in form of the pressure gradient balancing 
all external forces. In the following, we will therefore fo- 
cus the discussion on the case where the pervaporation 
vanishes before the diffusion coefficient diverges, 

(/)°^ < (j)"^. (73) 

For our treatment of q and D, this implies that the di- 
vergence of D is not reached, because q vanishes first. We 
thus expect the effect of drying out at least for very high 



densities. (Strictly speaking, this is true only for non-di- 
verging C{(f)).) In dense solutions, especially at the end of 
the pervaporation channel, we thus expect a flat concen- 
tration profile with a nearly vanishing pervaporation and 
a large but finite diffusion coefficient. 

We have to add a technical point here, as the diffusion 
coefficient depends not only on the osmotic compressibility 
but also on the phenomenological coefficient C{(j))- There 
are three different possibilities, which are all sketched in 
Fig. [51 If £(1) is non-zero, the diffusion coefficient diverges 
as ^ 1. This is the case which we have discusses so 
far. In the cases with vanishing phenomenological coeffi- 
cient C{1) = the diffusion coefficient remains finite or 
eventually vanishes. 



4 Prototype solutions: several phases 

After the treatment of smooth concentration profiles such 
as the one in Fig. [H we now continue with situations com- 
prising several distinct thermodynamic phases, see Figs.[9l 
and 1121 As more and more solute is transported from the 
reservoir into the channel the concentration rises and may 
lead to a phase transition. The position of such a phase 
interface and its movement in time is a convenient observ- 
able in experiments. The velocity of such an interface is 
governed both by the advective and the diffusive currents 
around it, and it thus presents a fingerprint of the mix- 
ture properties which have here been reduced to the two 
coefficients D{(j)) and q{4>). It does not provide as much in- 
formation on these coefficients as for example the full con- 
centration profile. The interesting question is how much 
information on the coefficients can be extracted from the 
velocity of an interface as a function of time. 

In addition to the model Eqs. (|T]) and ([2]) we now 
take the positions Xi{t) of the interfaces into account as 
unknown variables. This type of problem is known as a 
free-surface or Stefan problem. The interface movement is 
governed by Eq. (|39p . or in its dimensionless variant by 
Eq. (j66p . The typical situation with one moving interface 
is sketched in Fig. [SI 

4.1 Two phases: Moving phase interfaces 

To start with, we repeat the argument in Ref. [8J on the 
experiment of solution of AOT in water. It led to an expo- 
nential decrease of the interface velocity, or equivalently 
to an interface position of the form 

M = i_exp{-^^(t-to)} (74) 

with an effective nucleation time to. This result is based on 
several assumptions: first, that all incoming solute current 
is transferred into the interface movement. This means 
that the total solute current j{x,t) is a step- function, be- 
ing zero in the dense phase (left-hand side) and a non- 
zero constant in the dilute phase (right-hand side). Equa- 
tion (|74|l assumes further that the velocity profile vanishes 
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Fig. 6. Two different mechanisms to keep the density pro- 
file bounded from above: Panel (a) shows the influence of a 
diverging diffusion coefficient D'{(j)) — 1/(1 — </>) while the per- 
vaporation coefficient is constant (g' = 1, Pe = 10^). In this 
case, the velocity profile is linear for all times (not shown). 
In the two lower panels (b) the diffusion coefficient is constant 
and the pervaporation coefficient varies with tj) as q'{(j}) = l — cj) 
{D' — 1, Pe = 10^). This leads to a time-dependent velocity 
profile. 



everywhere left of the interface and that it is linear oth- 
erwise, v'^{x) = —Vwq{0)(^x — X{t)y With these two as- 
sumptions, the solute current which enters the channel at 
the reservoir (x = L) is 

j{L,t) = 4>Lv\L,t) = -v^q{Q)4>L[L - X{t)]. (75) 

The amount of solute coming from the reservoir thus de- 
pends on the position of the interface. As all solute current 
arrives at the phase interface, we find its movement from 




pervaporation 



refilling from outside 

7»(0) <7; 



Fig. 7. A sketch of the chemical potential 711, (<^) of water with 
the maximal concentration permitted by thermodynamics, (j)^'^. 
The pervaporation vanishes at the equilibrium value (jf'^ be- 
cause of the equilibrium '^^{(ff^) = 7^^"'. 



pefrheable 




, impermeable 
TP const 



impermeable 

Fig. 8. Three possible functions £(</>) for the phenomenological 
coefficient of inter-diffusion: They lead to qualitatively different 
behavior of the diffusion coefficient at </> — > 1. 



Eq. dMD to be 



m = - 



j{X{t),t) _ J^^q{0)cj>i 



A(j) 



A(j) 



{L-X{t)). 



(76) 



This differential equation leads to the solution ([7^ . 

The two assumptions are quite strong, since the tem- 
poral change of the volume fraction (f) does lead to a spatial 
change of solute current. Nevertheless, the solution (ffi)) 
has been successfully fitted to experimental data, at least 




Fig. 9. Sketch of the typical situation with one moving phase 
boundary, described in Sec. 14.11 
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for short times. For longer times, however, it over-esti- 
mated the velocity. It is this deviation which we like to 
explore in the following. 

For a more detailed analysis of the interface movement 
we have to develop a global view on the spatial solute 
profile. In general, whatever value the solute current has 
at the entrance of the channel (i. e. the concentration (jy^ 
times the velocity there) , it drops to zero at the end of the 
channel. We may therefore identify three major possible 
growth mechanisms due to the spatial change of the solute 
current along the channel: One possibility is the case we 
discussed around Eq. ([74|) , where the solute current is con- 
stant in both phases with a discontinuity at the interface. 
The current leads to the growth of the dense phase and to 
the movement of the interface. The two other possibilities 
are losses in the bulk of the dilute and the dense phases, 
leading to a temporal change of the concentration profiles 
there. 

Independent of the total solute current, the mixture ve- 
locity v'^ may be zero or not in the dense and dilute phases. 
In the following three subsections we discuss three dif- 
ferent combinations of vanishing/non- vanishing currents/ 
velocities left and right of the interface. 



4.1.1 Solute accumulation at the phase boundary: the 
spatial profile in the dilute phase 

We first investigate the spatial concentration profile which 
results from the assumptions taken above for Eq. (j74p . 
namely that all incoming solute is accumulated at the 
phase interface. The dense phase plays no role, i. e. v^{x) = 
0, j{x) = for X < X{t). The velocity of the interface 
then depends only on the concentration profile in the di- 
lute phase, 

M»-2^^ (77, 

In a system of coordinates moving together with the inter- 
face the evolution equation for the density reads (compare 
with Eq. (p]) ). 



5t0 = Xd^(t) - dxj. 



(78) 



A steady solution in this moving frame corresponds to a 
uniform current jo = —'^4>X in the reference frame. The 
uniformity of the current corresponds to the above men- 
tioned assumption that all solute entering the dilute phase 
leaves it at the interface. We then find the profile station- 
ary with respect to the moving frame as the solution of 



(79) 



As this equation is to be solved in the dilute phase, we may 
assume constant coefficients q{4)) k, q{Q) and D{4)) k, D{{)) 
which allows to find the spatial profile 



0(X(i)+y) =cxp 



2A4>- 

a 



G 



G 



(!) 



(80) 
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Fig. 10. Visualization of the solution (|80|) which is stationary 
in the frame moving with the interface. The parameters are 
arbitrary (Pe = 10^ X' = 0.05, = 0.3, 0*'^' = 0.8). 




13 , G{x) := / e^'di. 



with the shortcuts 

X 

Vwq 

(81) 

An example of this solution is shown in Fig. (TUl It exhibits 
the same two regimes as the dilute equation described 
above in Sec. l3.2l and in Fig.O There is a hump dominated 
by diffusion and a hyperbolic ramp. 

The solution from Eq. ((80|) . depicted in Fig. [TOl gives 
insight into the implications of the taken assumptions: It 
shows that the diffusive hump exceeds the binodal con- 
centration value (/)*•'''' . Evidently, there is a zone of higher 
concentration which moves together with the phase inter- 
face, and which continuously feeds the dense phase. As 
long as both the velocity and the total solute current van- 
ish at the interface, the diffusive hump must exceed the 
value (/)(''^ since there must be a negative current causing 
the interface to move to the right. The only contribution 
to the solute current is the one caused by the concentra- 
tion gradient, thus having to be positive. The assump- 
tions j{X^t) = and v^{X,t) — create a situation in 
which it is in principle possible that the peak which is 
clearly visible in Fig. [10] reaches a spinodal value, which 
is larger than 0'^'') but smaller than As the concen- 
tration reaches the spinodal, it will undergo local phase 
changes and exhibit local blobs of a metastable phase. 
Such blobs have occasionally been observed in the exper- 
iments using a KCl solution [T . 

The diffusive hump in Fig. [TU] appears pronounced due 
to the large Peclet number chosen in this example. The 
Peclet number renders the length scale a small, which re- 
duces the lateral extension of the hump and increases its 
maximum value. In experimental realizations, its height is 
limited by several factors, first by the overall Peclet num- 
ber, second by the spinodal concentrations which enforce 
a local phase transition instead of the smooth maximum, 
and third by the pervaporation in the dense phase. For a 
non- vanishing velocity v^{X,t) the necessity of having a 
positive slope dx(t> sd x = X{t) gradually vanishes, allow- 
ing the height of the maximum to decrease or to vanish 
completely. 
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We continue to discuss the solution ([80| by showing 
how the result (f74|) follows from the spatial form of the 
solute concentration. The hyperbolic ramp, which is the 
asymptotic solution for large x, here reads (in rescaled 
units) 

(t)(X(t)+y) (82) 

We may use the asymptotic form to calculate an approxi- 
mate velocity of the interface: Of course, this result corre- 
sponds to a stationary solution in a reference frame mov- 
ing at constant velocity. Therefore, the input concentra- 
tion at the end of the channel must be increased consis- 
tently in order to find this solution. Reversely, if the input 
concentration is held constant at (l){L) = (j)^, the current 
arriving at the phase interface decreases in time. As long 
as the diffusive hump does not reach the inlet of the chan- 
nel, this change of current is small and we may use 4>{L) 
as a slowly varying parameter. Setting 4'{L) — (f>L in the 
asymptotic solution (|82p then yields essentially the same 
expression for the evolution of the interface position as in 
equation ([74| . 



X'{t) = 1 - expl 



-it'-t',)}. (83) 



4.1.2 Solute accumulation at the phase boundary: the 
velocity field in the dense phase 

As a next approximation we relax the assumption of a 
vanishing velocity left of the interface. However, we con- 
tinue to keep the total current zero in the dense phase. 
Also the dilute phase is passive, such that all incoming 
solute is accumulated at the interface. 

As discussed above in section 15751 the velocity and den- 
sity profiles in dense phases can behave in two qualita- 
tively different ways: One corresponds to the vanishing 
of q{(t)) and the other to the divergence of D{(f)). We make 
the same distinction here. The first case, with vanish- 
ing 9 ((/)), is characterized by a flat velocity profile near the 
end of the channel. It is likely to lead to a density profile 
which has maximum density from the end of the channel 
up to the interface, where it reaches a small region with 
non-vanishing velocity. After a transient time, the den- 
sity profile in a small region could become stationary in a 
reference frame moving together with the interface. This 
solution supports a mixture velocity at the interface that 
is constant in time. The position of the interface is then 
given by 



m 

L 



= 1- 



Liytvq{0) 



)(l-exp{- 



Acj) 



The dilute phase is here characterized by a constant q and 
D. 

In the second case, in which the diverging diffusion 
coefficient leads to a bounded density profile, the veloc- 
ity field needs not to vanish except at the very end of 
the channel. For simplicity, we adopt a constant non-zero 



value for the pervaporation in the dense phase, (7(0) « 
g(l) < (7(0). This solution leads to a piecewise linear mix- 
ture velocity and to the interface position 



Mt) 

L 



{l-..^[-^-^it-t,)}). (85) 



9(0) -9(1) 



A(j) 



None of the solutions (jM)) and ([85)1 can explain the 
velocity decrease found in the AOT experiments. Instead, 
both velocities are larger than the initial proposition solu- 
tion (|74p . Note that we still kept the assumption j {x) = 
in the dense phase, which corresponds to a large diffusion 
coefficient. It must therefore be concluded that in the case 
of AOT either the diffusion coefficient is not sufficiently 
large, or the velocity decrease stems from the dilute phase. 



4.1.3 Solute accumulation also in the dilute phase 

In both preceding subsections we have kept the assump- 
tion that all incoming solute contributes to the movement 
of the interface. We now continue with the more general 
case where there may be a temporal change of the con- 
centration profile such that there is a gradient of solute 
current in the dilute phase. The dense phase is still as- 
sumed to be passive (no pervaporation). 

This more general situation cannot be described ana- 
lytically but has to be treated numerically. A numerical 
test of the situation with a dense phase growing into a di- 
lute one exhibits the limitations of the assumption taken 
above in Sec. I4.1.D Figure [TT] shows the full numerical so- 
lution of the model. The velocity in the dense part has 
been kept zero [q — 0), leading to a flat density profile 
there. The dilute part of the solution agrees qualitatively 
with the graph given in Fig. 1101 There are, however, devi- 
ations. The plot shows the solutions of Eq. ([80|l with the 
true position and the true velocity of the interface taken 
from the numerical solution. Panel b quantifies the differ- 
ence between the approximation from Sec. 14.1.1! and the 
true solution. It shows the position of the interface as a 
function of time. A fit with the exponential function (j74p 
with the known value for (jj^iy^q/ A(j) and to as the fitting 
parameter proves that the approximation does not work 
in this example. A second fit for both parameters indi- 
cates that the resulting function is still exponential, but 
with a different timescale than assumed in the approxima- 
tion ([74)1 . In both fittings only the data in the gray-shaded 
area has been used. The fit quantifies that around 60% of 
the incoming solute is accumulated at the interface and 
leads to the growth of the dense phase. The rest remains 
in the dilute phase. Panel fTTb visualizes this ratio directly. 
Evidently, the current has not the form of a step function, 
but takes on different values at the interface and at the 
channel entry. This clearly falsifies the assumption of a 
uniform current in the dilute phase. 

Upon varying the parameters of the example, we found 
that the concentration 0^''^ right of the interface has the 
largest infiuence on the deposition of current. For = 
0.05 nearly all of the incoming current was found to ac- 
cumulate at the phase boundary — despite a well devel- 
oped diffusive hump similar to the one in Fig. Illb . and 
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Fig. 11. Numerically determined density profile cf>{x' ,t') con- 
taining a single moving phase interface. The binodal density 
values at the interface are 0'*-* = 0.8 and = 0.3; at the 
channel entry (f)L ~ 10~^; Pe — 10'^. Panel (a) shows the den- 
sity profile, panel (b) the position of the interface as a function 
of time, and panel (c) the evolution of the spatial current pro- 
file. 



despite the fact that the solute current still exhibited a 
pronounced hump. The current apparently overflows the 
diffusive hump completely, if t) has already reached 
the value 0^'') before the hump starts. We may quantify 
this by stating that approximation (I74p is valid if 



< 



Pe 



L-X 



(86) 



Here, we assumed the width of the diffusive hump to scale 
as L/V Pe , see the definition of a in ([STjl . and also Ref. T. 



4.2 Several phases: Phase slab thicknesses in 
stationary solutions 

We now return to the introductory example of Sec. 11.2! 
with four phases, one of them dilute. This situation is 
sketched in Fig. [121 The aim here is to present an example 
how the thicknesses of phase slabs are governed by the 
thermodynamic and dynamic properties of all the phases. 

For the solution of AOT in water, a remarkable ob- 
servation has been reported in Ref. [8]: First, the same 
phases as in equilibrium have been found in the evapora- 
tion channel. They were identified by their polarizability 
and birefringence, showing hexagonal, cubic and lamellar 
internal structures for the dense phases, while the dilute 
phase was isotropic. The spatial pattern showed that the 
two enclosed slab thicknesses (of the cubic and the lamel- 
lar phases) were of comparable extension. This result is 
remarkable because these two phases occupy a very dif- 
ferent range of 0- values in the equilibrium phase diagram. 
There, the lamellar phase appears much broader than the 
cubic phase, see the supplementary material of Ref. [8]. 
We conclude that there is apparently a much larger con- 
centration gradient in the lamellar phase than in the cu- 
bic phase. This larger concentration gradient was stable 
in the experiment, even when all dense phases were con- 
tinuously fed with more incoming solute. After a short 
transient time, the interfaces bounding the cubic and the 
lamellar phase moved together at constant distance. 

We now show that our model possesses qualitatively 
the same solutions. We also observe stable thicknesses of 
comparable size for the lamellar and the cubic phases, if 
we choose the chemical potential and the phenomenolog- 
ical coefficients well. The strategy here is to adopt simple 
model functions for 7u,((/>) and £(</>). They are then trans- 
ferred into diffusion coefficient D((f>) and pervaporation 
coefficient q{(l)) using Eqs. ([61]) and (j62|) which arc then 
used in the stationary model equations. 

We here adopt the most simple models for the chem- 
ical potential and for the phenomenological coefficients. 
The functions used below are displayed in Fig. [131 The 
values of (j) are taken from the equilibrium phase diagram 
of AOT (see supplementary material of Ref. [8, ) . The re- 
gions shaded in gray indicate phase coexistence, where 
the chemical potentials are constant. Within the most 
dense phase we use a logarithmic behavior of the form 
7^^,(0) oc ln(l — (j)), whereas in the other three phases we 
linearly connect some coexistence values which have been 
chosen as model parameters. The phenomenological coef- 
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Fig. 12. Sketch of the situation with several phase slabs ex- 
hibiting characteristic extents, as described in Sec. 14.21 
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ficient has been chosen constant in each phase except in 
the dilute phase where a linear function of is required 
(see Sec. [O^). The coefficients and D'{(p) follow 

from equations ((6T|) and ([3T|l , see Fig. [TH 

Figure (TS] shows a typical profile of the volume frac- 
tion of solute, 4>{x'), as it results from the model functions 
given in Figs. [13] and [Ml The lamellar and cubic phases 
are indeed found to have approximately the same extent. 
It is mainly the different slope of the chemical potential 
in Fig. [T3k which has been adjusted manually to achieve 
the comparable extent of the two phases. The physical 
content of these slopes are the very different osmotic com- 
pressibilities of these two phases. Of course, also the phe- 
nomenological coefficient plays a role. We do not claim 
that the presented chemical potentials are the only ones 
which may lead to comparable phase extensions. However, 
it is evident in this very example that the osmotic com- 
pressibilities of the chemical potential directly influence 
the extension of the phase slabs. The profile in Fig.[T5lhas 
been calculated as the stationary solution of the model 
equations ([SG]) and (jST]) . The use of the stationary equa- 
tions will be explained in the following. 

4.2.1 Stationary and quasi-stationary solutions 

The full evolution of the phase thicknesses in time is sub- 
ject to many parameters including not only the model 
parameters for D{(j)) and q^cj)) and the incoming volume 
fraction and the Peclet number. It depends also on the 
initial state, or in other words, on the precise conditions 
of phase creation at the end of the channel. We do not 
want to treat the issues of phase creation here, the more 
as a reasonably fast numerical method for this parameter 
regime is not at hand. Instead, we employ a quasi-station- 
ary treatment, which is focused on the main features of 
the phase thicknesses. The idea is to assume a very small 
input current and to assume at any time an equilibrated 
stationary profile. Such a stationary profile is shown in 
Fig. [15] for the model functions of Fig. [13] Note the slab 
thicknesses of the lamellar and the cubic phases which are 
of the same extent. This behavior has been achieved by 
choosing the slopes of the chemical potential in Fig. [TS] 
which is 20 times larger in the cubic phase than in the 
lamellar one. Of course, the variation of the chemical po- 
tential also has its influence on the evaporation. We thus 
cannot clearly separate between thickness due to a varia- 
tion of diffusion and due to modifled pervaporation. This 
ambiguity is a systematic issue in the interpretation of the 
experimental results, and it is reflected in our model by 
the two coefficients D and q. 

In order to obtain the ensemble of all boundary po- 
sitions, which is given in Fig. 1161 we parametrized the 
solutions by the total amount of solute / 4){x')dx' . In a 
quasi-stationary setting with extremely small input cur- 
rent, this amount may be understood as time multiplied by 
the input current. A numerical test of the full time-depen- 
dent equations revealed similar lines as a function of time 
(not shown). It can be seen that the interior slab thick- 
nesses relax to constant values, and that only the last, the 
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Fig. 13. Simple model functions for the rescaled chemical po- 
tential difference of water across the PDMS layer and for the 
phenomenological coefficient, both as functions of the volume 
fraction 4>- Note the different slopes of the chemical potential 
in the lamellar and the cubic phases which lead to the same 
slab thickness in Fig. [I5l 
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Fig. 14. The rescaled coefficients for evaporation q [cj>) and for 
diffusion D'{(j)) resulting from the model in Fig. 1141 

hexagonal phase continues to grow. This behavior stems 
from the vanishing pervaporation in the hexagonal phase 
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Fig. 15. An example solution of the stationary equations in 
which the thickness of phase slabs is controlled by the slope of 
chemical potential in Fig. [TSl Note the equal thickness of the 
lamellar and the cubic phases. 
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Fig. 16. The development of phase thicknesses for (a) the 
model given in Fig. 1131 and in (b) for a chemical potential 
in the cubic phase having the same slope as the one in the 
lamellar phase. 



which leads at the interfaces to a mixture velocity con- 
stant in time. If the velocity does not vanish in the dense 
phase, then the interior thicknesses do not truly become 
stationary. 

In order to demonstrate the effect of the slopes of 
the chemical potential on the slab thicknesses, Fig. [TBb 
presents the result for a modified model chemical poten- 
tial. We took the curve in Fig. [T3k and changed the slope 
in the cubic phase to the slope in the lamellar phase. 
This situation implies the same osmotic compressibilities 



in both phases, and also a comparable spatial concentra- 
tion gradient. As a result, the cubic phase has shrunk and 
is hardly visible in Fig. 116b . 

The numerical solutions in Figs. fTSl and fTHl have been 
created by a Runge-Kutta algorithm solving the station- 
ary equations ([T]) and ([2]). On the analytical side, they may 
be solved formally by integration over (p. The stationary 
equations can be rewritten such that 

— [v") = v^q[(t>) and (87) 

ax cp 

dx _ D{(p) 

d(j} v^<f> 

The formal solution of the problem then reads 



u°(a;,xo) 



K(-o))' 



2v^ J 



(89) 
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(90) 



This expression describes in principle any stationary 
solution, but in the unusual way x{(p). This formulation is 
especially convenient for the determination of phase slab 
thicknesses. For the first phase at the end of the evapora- 
tion channel, the expression simplifies due to the bound- 
ary condition u°(0) = 0. In case of the phase interface 
X = X(t) the lower bound of the integrals is fixed to the 
binodal value 0^^' of the phase transition. We then find 



0(0,4) 



(91) 

L 
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4.2.2 Phase nucleation 

The treatment of the phase thickness in the present sec- 
tion resembles to some extent the treatment of nucleation 
as it has been presented by Evans, Poon, and Gates |18j . 
Both include the assumption of local equilibrium, leading 
to a diffusion-type equation, while the total system is out 
of equilibrium. Both allow the occurrence of metastable 
phases which are absent in equilibrium. Here, these phases 
are the dense ones growing because of the non-zero solvent 
flow; in Ref. [TS] it is a freshly nucleated metastable phase 
which grows because the whole system is quenched in tem- 
perature. 

We have not included the details of phase nucleation 
near the end of the channel, as it would require a more 
advanced numerical treatment. The general framework, 
however, should be applicable to the nucleation in a sim- 
ilar manner as in Ref. 1181. 
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5 Summary 

We have derived a model for the behavior of binary solute- 
solvent mixtures in recently developed microevaporator 
devices. The final target of the model is to allow an inter- 
pretation of experimentally observed data in terms of well- 
defined thermodynamic and dynamic properties of the 
mixture which are measurable also by different means. As 
a first step towards this goal we have developed the model 
equations ([!]) and ([2]) which relate the phase behavior in 
the microevaporator to an inter-diffusion coefficient D and 
another coefficient q quantifying the pervaporation of so- 
lute. We have shown in Eqs. (|3T|) and (|46|) how these two 
coefficients depend on the chemical potentials in the mix- 
ture and on the phenomenological coefficient for inter-dif- 
fusion. This relation has then been extended to make also 
the connection to other quantities such as the coefficients 
of sedimentation, permeability and osmotic compressibil- 
ity. 

The presented model is focused on the transport equa- 
tions of solute and solvent mass in binary mixtures, using 
the hypothesis of local thermodynamic equilibrium. We 
provide in detail all the necessary assumptions to reduce 
the description to two variables, one of which is the volume 
fraction of solute, the other one is the volume-averaged 
mixture velocity. We further assume incompressibility of 
the mixture and reduce the spatial dimensions to the one 
parallel to the evaporation channel. 

We solved the model equations in a number of typical 
limit cases which still allow analytical approximations and 
compared these with full numerical solutions. In particular 
we found good agreement in the case of dilute solutions, 
where the diffusion and pervaporation coefficients may be 
treated as constants. In dense binary mixtures, where this 
is no longer the case, we were able to show that there are 
two qualitatively different mechanisms which both lead to 
a saturation of the concentration at its maximal possible 
value. One of them is due to vanishing evaporation at an 
equilibrium concentration. The other case comes from the 
singular behavior of the chemical potential of water in this 
limit. 

We have further analyzed the growth velocity of a 
dense phase into a dilute one and have found several ex- 
pressions which allow to extract quite some information on 
the properties of the mixture only from tracking a single 
interface between two phases. Further, we have analyzed 
the resulting thicknesses of phase slabs in the microevapo- 
rator in a stationary setting. In particular, the experimen- 
tally observed case of equal phase thicknesses has been 
regarded. We were able to explained this remarkable fact 
as the result of very different osmotic compressibilities in 
the two phases. 

In total, the results presented here show that the in- 
terpretation of the experimental results is possible but far 
from trivial. In all cases, there is some ambiguity in the 
interpretation of the observed behavior: As the physical 
properties of the mixture are summarized in two coeffi- 
cients, also the reasons for the flattening at high solute 
concentrations are twofold. In a similar manner, the thick- 
nesses of the phase slabs are not only governed by their 



thermodynamic properties, but also by the phenomeno- 
logical coefficient containing dynamic information. 

It is precisely this twofold information contained in 
the microevaporator data as well as in the two coeffi- 
cients D and q which presents the core of the present 
work. We are convinced that the presented model estab- 
lishes a good first compromise between the complexity of 
the phenomenon being out of equilibrium and the need 
for simplicity in the equations. Our model will help on 
the experimental plan to identify sets of measurements 
which allow to lift the ambiguity in the coefficients. It is 
the detailed derivation of the model and the numerical 
treatment of the equations which allows to improve the 
interpretation of the experimentally measured quantities 
and to extract from them the chemical potentials and the 
transport coefficients as full functions of the solute con- 
centration. 
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the viscosity of water, rj — 10 '^kg/(ms). This leaves us 
with a pressure gradient of 

^^r,a(w/h)^^3.9— (93) 

with the shape-factor a(lOO) ~ 1200 taken from Ref. [19] . 
As a further approximation (B) we take the chemical po- 
tentials in their dilute limit, where 7s becomes logarithmic 
in Cs, yielding 

(94) 

dCa nis Cs 

We also anticipate (C) the "hyperbolic ramp" solution at 
the entrance of the channel (see Sec. 13.21 and Ref. [7] for 
an experimental evidence), Cs{x,i) ~ Cs{L)L/x. This ap- 
proximation leaves us with the concentration gradient of 
the order |Vcs| ~ Cs{L)/L. Putting all pieces together, the 
ratio of the two terms in Eq. ((29} is 

{i's-u.u,)\VP\ . Lm. Pa 

dcs I I 

If we say that the solute is ten times lighter than the sol- 
vent water (approximation D), the difference of specific 
volumes will be dominated by that factor. Together with 
a molecular weight of a few hundred atom units, we end 
up with the factor in Eq. (|95|l to be less than 10~^, which 
justifies the approximation in Sec. 12.1.41 

15. W.H. Press and S. A. Teukolsky and W. T. Vetterling 
and B. P. Flannery, Numerical Recipes in C (Cambridge 
University Press, Cambridge 1992). 

16. G. C. RandaU, P. S. Doyle, Proc. Natl. Acad. Sci. USA. 
102, (2005) 10813. 

17. The pervaporation may cease if the spheres exhibit some 
compressibility which enters the chemical potential. 

18. R. M. L. Evans, W. C. K. Poon and M. E. Cates, Europhys. 
Lett. 38, (1997) 595. 

19. N. A. Mortensen, F. Okkels, and H. Bruus, Phys. Re- 
view E 71, (2005) 057301. 



